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ABSTRACT 

The fragmentation of star-forming interstellar clouds, and the resulting stellar initial mass function 
(IMF), is strongly affected by the temperature structure of the collapsing gas. Since radiation feed- 
back from embedded stars can modify this as collapse proceeds, feedback plays an important role in 
determining the IMF. However, the effects and importance of radiative heating are likely to depend 
strongly on the surface density of the collapsing clouds, which determines both their effectiveness at 
trapping radiation and the accretion luminosities of the stars forming within them. In this paper we 
report a suite of adaptive mesh refinement radiation-hydrodynamic simulations using the ORION code 
in which we isolate the effect of column density on fragmentation by following the collapse of clouds 
of varying column density while holding the mass, initial density and velocity structure, and initial 
virial ratio fixed. We find that radiation does not significantly modify the overall star formation rate 
or efhcicncy, but that it suppresses fragmentation more and more as cloud surface densities increase 
from those typical of low mass star-forming regions like Taurus, through the typical surface density 
of massive star-forming clouds in the Galaxy, up to conditions found only in super star clusters. In 
regions of low surface density, fragmentation during collapse leads to the formation of small clusters 
rather than individual massive star systems, greatly reducing the fraction of the stellar population 
with masses > 10 Mq. Our simulations have important implications for the formation of massive 
stars and the universality of the IMF. 

Subject headings: ISM: clouds — radiative transfer — stars: formation — stars: luminosity function, 
mass function — turbulence 



1. INTRODUCTION 

The stellar initial mass function (IMF) is determined 
by the fragmentation of gas clouds into progressively 
smaller pieces as they collapse, and the characteristic 
fragment mass scale is thought to be proportional to 
the Jeans mass. However, in an infinite, isothermal gas 
cloud the Jeans analysis does not pick out a unique mass 
scale for fragmentation. Since the Jeans mass varies as 
Mj (X p~^/^T'^/^, where p and T are the gas density 
and temperature, in an isothermal cloud Mj reaches ar- 
bitrarily small values as the cloud collapses and p in- 
creases, allowing fragmentation to proceed to arbitrar- 
ily small scales. In terms of numerical simulations, this 
lack of a natural fragmentation scale in isothermal clouds 
is reflected in the fact that the amount o f fragmenta- 
tion i s ultimately resolution-dependent (e.g. Martel et al. 
20061, and that an isothermal simulation in a periodic 
box can always be rescaled so as to give the fragments 
that form an arbitrary mass, while maintaining a fixed 
virial parameter, Mach number, and number of initial 
Jeans masses (e.g. [Offner et al.]|2008a[ ) . 

For this reason, any explanation of the fragmentation 
scale of molecular clouds and thus the stellar IMF re- 
quires a deviation from uniform, isothermal flow. How- 
ever, the nature of this deviation remains controversial. 
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In tu r bulent fragmentation mo dels fPadoan & Nordlund 
2002| JHennebelle fc Chabrier| |2008 , 2009), the shape of 



the LMb' IS determined by the properties of turbulence, 
but the overall mass scale is set by computing the Jeans 
mass at a density proportional to the mean densitjr] ei- 
ther in an entire giant molecular cloud or in some smaller, 
fiducial star-forming gas clump. While there is signifi- 
cant support from nu merical sim ulations that a process 
of this sort operates ( [Padoan et al. 2007), and observa- 
tions in dicate that the Jeans scale is imprinted in CO 
clumps (Blitz & Williams 1997), this model has two sig- 



nificant gaps, l^'irst, since molecular clouds have com- 
plex structures that span a large range of densities, it is 
not obvious how to define the star-forming region over 
which the mean density should be computed. Second, 
once a gravitationally-bound collapsing object is formed 
in these models, it is unclear why it should not fragment 
even further as it collapses, since the bound object now 
defines a new cloud with a higher mean density, and thus 
a smaller fragmentation scale. Indeed, purely hydrody- 
namic simulat ions of isolated ma ssive cores find exactly 
this behavior ( |Dobbs et ar] |2005'). 

In contrast, in non-isothermal fragmentation models 
the characteristic mass scale is introduced via a small 
deviation from isothermality that occurs at some den- 
sity, which then sets the density and temperature that 
enter into the Jeans mass. The necessary kink in the 
equation of state may arise in several ways. It can come 

^ In models of this sort the relevant density is usually taken to be 
the density time s the square of the Mach number, which Krumholz] 
I&: McKee| | |2005"l l point out is equal to the Jeans mass computed at 
a pressure equal to the mean ram pressure. 
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from from a transition between optically th in and op 
tically thick conditions as gas collapses ( e.g. Masunaga 
eTadp^QSl [Masunaga fc Inutsuka||2000l [Jjonneil et aT 
20061) ; the opacity limit f or frag mentation, first proposed 
5y"Eow & Lynden-Bell (19761, is one example of such 
a transition, albeit at a mass scale of 0.004 Mq (Whit- 



worth et al. 2007), too low to be relevant for the bulk 
of stars. A second possible origin for non-isothermality 
is the density-dependent interaction between mole cular 
cooling, cosmic ray heating, and dust-gas coupling ( Lar 



son 



2005; Elmcgrccnet al. 2008). That non-isothermality 
ot this soft can aftect frag mentation also has su pport 
from numerical simulations ( [Jappsen et al. 2005|, but 
these models too face difficulties. I'heir procedure for 



procedr 

computing the characteristic density and temperature 
relies on an effective equation of state based on aver- 
age rates of radiative heating and cooling. This ignores 
the large spatial and temporal variations in the radiation 
field in star-forming clouds. 

These difficulties have led to a renewed focus on the 
potential importance of another mechanism for setting 
a characteristic mass scale: radiation feedback. Young 
stars radiate prodigiously, due to accretion at low masses 
and via Kelvin-Helmholtz contraction and nuclear burn- 
ing at higher masses, and these effects can heat the 
gas around them , rais ing the Jeans mass. Analyt ically, 
Krumholzl ( |2006[ ) and [Krumholz fc McKee| ( |2008[ here- 
after KM08) show that stellar feedback should strongly 
suppress fragmentation by raising the temperature in 
star-forming clouds. Numerical simulations of high-mass 
star (Krumholz ct al. 2007a, hereafter KKM07) and low- 
mass ( jUftner et al.||2009b| ) star formation confirm this 
conclusion. Under some circumstances feedback effects 
completely swamp the subtle changes in the equation of 
state on which the non-isothermal fragmentation models 
rely. 

While the importance of radiation feedback has been 
recognized, its relative importance in different star- 
forming environmen ts is on ly starting to be considered. 
Analytic models by KM08 suggest that effectiveness of 
feedback will depend on the column density of the star- 
forming cloud, which determines its ability to trap pro- 
tostellar radiation. The effect is strong only when S > 1 
g cm~^. Observed star- forming clouds have surface den- 
sities ranging from S 0.1 g cm~^ in diffuse c louds such 
as Perseus and Ophiuchus (Evans et al. 2009) to S ~ 1 
g cm~^ in typical Galactic regions of massive star for- 
matio n ( Shirley et al.||2003 Faundez et al.|2004 Fontani 



et al. 2005) to E ^ 10 g cm ^ o r more i n extragalactic 



super star clust ers ([Turner et a l. 2000; McCrady & Gra- 
cK " 



ham|20 07'), and |Mckee fc 'l'an[ (|2'002n2003| hrst pointed 
out that massive stars seems to form only at the high sur- 
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face density end of this distribution. However, numerical 
simulations thus far have not explored this parameter 
space systematically. [KKM07[ find that feedback sup- 
presses the fragmentation of niassive pro tostellar cores 



with E ~ 1 g cm-2, while Offner et al. (2009b) find a 



similar but significantly weaker suppression of fragmen- 
tation at E ~ 0.1 g cm~^. 

The only comparisons of fragmentation with radiation 
in clouds of varyi ng surface density performed t hus far 



are those of Bate (2009) and Urban et al. (2010). Both 



well below [KM08[ s analytically-predicted threshold for 
fragmentation, and well below the observed column den- 
sity in the typical region of star cluster formation in the 
Galaxy. Furthermore, neither include radiative trans- 
fer in a way that is appropriate to study suppression of 
fragmentation in dense regions. In Bate's simulations, 
gas can radiate only up to the point where it is captured 
by a sink particle. The radiation it releases when it ac- 



cretes o nto the stellar surfac e is ne glected. [Offner et al 
([2009b') find that this leads |Bate[ to underestimate the 



of these authors consider surface densities ^ 0.4 g cm 



energy budget available for heating the gas by a factor 
of 20. Urban et al.'s simulations include radiation from 
stars but not radiation produc ed by either co mpression or 
viscous dissipation, although Offner et al. s results sug- 
gest that this approximation is not bad once stars be- 
gin heating the gas. However, Urban et al. also rely on 
a pre-computed spherically-symmetric profile that does 
not account for deviations from spherical symmetry in 
the surrounding density field. This is reasonable for the 
low optical depth clouds and relatively large length scales 
they consider, but it is questionable for dense cores and 
on small spatial scales, where non-symmetric shielding 
effects can be important - for example accretion disks 
are typically colder than the gas above or below them, 
due to their large optical depths, and Urban et al.'s ap- 
proach would not capture this effect. 

The goal of this paper is to fill that gap by studying 
molecular cloud fragmentation and the stellar mass func- 
tion in different star-forming environments while self- 
consistently taking into account the effects of stellar ra- 
diation feedback. We perform a controlled experiment 
by running a series of simulations of collapsing, turbu- 
lent gas clouds in which we hold fixed the initial cloud 
mass, temperature, virial ratio, and turbulent velocity 
field, while varying the cloud surface density. To ensure 
that the results are robust, we are careful to hold nu- 
merical aspects of the calculations fixed as well. Each 
simulation uses the same numerical method, the same 
criteria for refinement, and the same maximum resolu- 
tion, so that any differences in outcome should be solely 
the result of radiative effects. In Section [2] we describe 
the numerical method and initial conditions we use for 
these simulations. In Section|3]we report the results, and 
study how fragmentation varies with initial surface den- 
sity. Finally, in Section |4] we discuss the implications of 
our results for massive star formation and for the IMF 
more generally, and we summarize in Section [sj 

2. NUMERICAL METHOD AND INITIAL CONDITIONS 

2.1. Equations and Solution Algorithms 

Our nume rical method is identical to that of lKrumholzl 
et al. (j2009 ) , and we give an extensive description in the 
Supporting Online Material of that paper, so we only 
summarize our methods briefly here. Our simulations use 
the parallel adaptive mesh refinement (AMR) radiation- 
hydrodynamics code ORION. In ORION, the gas plus ra- 
diation fluid is represented at every grid point by a vector 
of conserved quantities (p, pv, pe, E), where p is the den- 
sity, V is the velocity, e is the specific non-gravitational 
energy (including kinetic and thermal), and E is the ra- 
diation energy density. The computational domain also 
contains an arbitrary number of point mass "star" par- 
ticles, each of which is characterized by a position a;, 
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a mass M, a momentum p, and a luminosity L. The 
code updates these quantities by solving the equations 
of radiation- hydrodyna mics plus gravity in the conserva- 
tive, mixed-frame form (Mihalas & Klein 1982), retaining 
terms to order v/c accuracy, and using the tlux- limited 
diffusion approximation to represent the radiation flux 
( Krumholz et al.||2007b[ ). The equations for the gas are 



d 



— p = - V • (pv) - M^W{x - X,) 



d_ 
dt 



{pv) = -V • {pvv) - VP - p\/(t> - WE 



(1) 



(2) 



-V • [{pe + P)v]~ pv-V(j)- Kopp(47rB - cE) 
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vE 



Xi), 



(4) 



where the summations run over all star particles present, 
Mi, and £i are the rates at which mass, momentum, 
and energy are transferred from gas to star particles, and 
W{x — a^) is a weighting function that distributes the 
transfer over a kernel 4 cells in radius. We calculate these 



using the Krumholz et al. ( 2004 1 sink particle algorithm, 
which we summarize below, i'he corresponding evolution 
equations for the star particles are 



-M,=Mi 



dt 



dt 
d 

dt^^ 



Pi. 



In these equations, the gravitational potential 

by 



= -47rG 



YMiS{x- Xi) 



The pressure P is given by 
P = 



(5) 
(6) 

(7) 
is given 

(8) 
(9) 



where Tg is the gas temperature, p, — 2.33 is the mean 
molecular weight for molecular gas of Solar composition, 
and 7 is the ratio of specific heats. We adopt 7 — 5/3, ap- 
propriate for gas too cool for hydrogen to be rotationally 
excited, but this choice is essentially irrelevant because 
Tg is set almost purely by radiative effects. The remain- 
ing quantities are the comoving frame specific Planck- 
and Rosseland-mean opacities kqr and kqp, the Planck 
function B = cauT^ / (47r) , and the flux limiter A and 



Eddington f actor R2 , computed using the [Levermore fc:| 
Pomraning ( 1981 ) approximation: 



A = — I cothi?- — 
R 



R V 

KORPE 

R2^X + X^R^ 



R = 



(10) 

(11) 
(12) 



We obtain the dust opacities K gp and kqb. f rom a 



piecewise- linear fit to the m odels of Pollack et al.| (|1994[) 
see Krumholz et al. ( 2009 1 for the exact tunctiona 



lor 



It is worth noting that the Rosseland opacity we use in- 
cludes absorption but not scattering effects, and as a 
result is likely something of an underestimate. Using a 
higher opacity would likely enhance the radiation effect 
we describe below, as sugge sted by recent s tatic radiative 
transfer calculations ( jPunham et al.p010[ ). 

We solve these equations in tour steps. Firstjthe hy- 
drodynamics module updates equations ([l} - (|3]) using 
all the terms on the right hand side except those involv- 
ing star particles or radiation. The update is based on a 
conservative Godunov scheme with an approximate Rie- 
mann solver, and is second- order a ccurate in time and 
space ( [Truelove et' al. 1998; |Klein| [r999). Second, the 
gravity module solves the Poisson equation to up- 
date the gravi tational potential using a rnultigr i d iter 



2002| . Third , the radiation module updates the right 
id si 



ation scheme (Truelove et al. 1998 Klein 1999 Fisher 



hancf sides of equations ([2| - ([4 ) for the terms in volving 



radiation. The module uses the Krumholz et al. (2007b I 
operator splitting method, in wnich the dominant terms 
describing radiation-gas energy exchange and radiation 
diffusion are updated using a fully i mplicit method based 
on ps eudo-transient continuation (Shestakov & Offner 
2008 1, and then the sub-dominant work and advection 



terms are handled explicitly. The fourth step is the star 
particle module. This portion of the code updates equa- 
tions ([1]) - ([7| using the gas-particle exchange terms. 
The code computes Mi by fitting the flow in the vicinity 
of each particle to a Bondi-Hoyle flow, and then set 
and £i by requiring that, in the frame comoving with the 
particle, accretion not alter the radial velocity, angular 
momentum, or temperature of the gas Krumholz et al. 
(20041. The method then updates the luminosity and 
the internal state of each star based on a simple one- 
zone protostellar evolution model. Details of the model 



are given in the Appendices of Offner et al. 
All of these modules operate with an AM 



(|2009b|). 

framework 



( Berger fc Ohger|1984| [Berger fc Collela|1989[ [Befl et al 



994|) in which we discretize the computational domain 



onto a series of levels I = 0,1,2, . . . , L. The coarsest level 
is Z = 0, which covers the entire computational domain. 
Subsequent levels cover sub-regions of the computational 
domain which are described by the union of rectangu- 
lar grids. Levels are nested such that any grid on level 
/ > must be fully contained within one or more grids 
of level I — 1. The cell size on level I is Ax;, and the 
spacings on levels I > are related to that on level by 
Ax I = Axo/^K The process of advancing the computa- 
tion through these levels is recursive. We first advance 
the grids on level by a time Ato, and then we advance 
the grids on level 1 by two timesteps of size Ati = Ato/2. 
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However each of these advances is followed by two ad- 
vances on level 2, and so forth, such each level advance 
involves 2' advances of timestep Aio/2' on level At 
the end of each advance of level I > through 2 steps, 
we perform a synchronization procedures between levels / 
and l — l to ensure conservation of mass, momentum, and 
energy across the interface between the two levels. We 
set the overall timestep Atg by computing the Courant 
condition (including a contributi on to the effective signal 



speed from radiation pressure - Krumholz et al. 2007b I 



separately on each level at the beginning of each coarse 
timestep. We then set AtQ — min(2'At'), ensuring that 
each level obeys the Courant condition for its advances. 

2.2. Refinement and Boundary Conditions 

The ORION AMR framework automatically adds and 
removes higher resolution grids throughout a simulation. 
We determine when higher resolution grids are required 
based on the following criteria: 

1. Any cell with a density greater than half the initial 
density at the edge of the cloud (see below) must 
be refined to at least level 1. 

2. We refine any cell within whose distance d to the 
nearest star particles is less than 16 Aa;;. This en- 
sures that regions around stars are always refined 
to the maximum level. 

3. We refine any cell where the density exceeds the 
Jeans density, given by 



pj = r 



TTCZ 



(13) 



where Cj is the sound speed and we use J = 1/16 



This avoids artificial fragmentation (Truelove et al. 
19971). 



4. We refine any cell where the gradient in the radia- 
tion energy density satisfies 



|V£;| > 0.15 



E 



(14) 



This ensures that we adequately resolve gradients 
in the radiation energy density. 

These conditions are applied recursively to every level up 
to some pre-specified maximum level L. In practice, the 
fourth condition is usually the most stringent. 

For all our calculations we use symmetry boundary 
conditions on the hydrodynamics, but we use adaptiv- 
ity to remove the boundary far enough from the cloud 
we are simulating so that no part of it ever approaches 
the boundary. The gravity module uses Dirichlet bound- 
ary conditions on the potential, with the potential on 
the boundary set equal to the value obtained from an 
octopole expansion of the density distribution inside the 
computational domain. The radiation module uses Mar- 
shals; boundary conditions, meaning that radiation en- 
ergy generated within the domain is free to escape. The 
incoming radiation flux is set to that appropriate for a 
blackbody radiation field at a temperature of 20 K. 



2.3. Initial Conditions and Simulation Setup 

The initi al setup f or all of our simulations is similar 
to those of |KKM07[ In all cases the setup consists of 
an initially spherical cloud of mass M = 100 Mq , initial 
mean surface density S, and radius R = ^ M / (ttE) in 
the center of a cubical computational domain. Observed 
dense star-forming clumps of molecular gas have roughly 
powerlaw density str uctures p oc r~''p with fc^ ^ 1.5 and 



considerable scatter (ICaselli & Myers|1995t|Beuther et al 

mm. 



2002 



20051 |2006l pueller et ai.lPO()j^l Itjridbaran et a 



20051 , so, folWng lMcKee fc Tan| ( |2003p , we adopt a 



powerlaw density structure with kp = 1.5 for our initial 
conditions. We give our clouds an initial turbulent ve- 
locity field chosen to put them in approximate balance 
between gravity and turbulent ram pressure. The veloc- 
ity dispersion is 



GM 



2{kp-\)R' 



(15) 



and the initial velocity field we use is identical to that 
of run lOOA of KKM0 7, which is a Gaussian-random 
field with a power spectral density P(k) cx fc~^ 



Fi- 



nally, we give the gas in the cloud an initial temperature 
Tg = 20 K, and we set the initial radiation energy density 
throughout the computational domain to E — 1.2 x 10~^ 
erg cm^"^, the value for a blackbody radiation field at a 
temperature Tr — 20 K. 

Outside the cloud we place a hot, diffuse ambient 
medium with a density pa — Podgc/100, where 



Pcdgc 



47r 



M 

i?3 



(16) 



is the density at the cloud edge. The ambient medium 
has a temperature Ta = lOOTg, ensuring that it is in 
thermal pressure balance with the cloud. To ensure that 
the ambient medium does not cool, radiatively heat the 
cloud, or interfere with radiation escaping from the cloud, 
we set the opacity of the ambient medium to a numeri- 
cally small value. 

We simulate three different clouds, chosen with values 
of E to be representative of three different types of star- 
forming environment as discussed in Section [T] The first 
run has S = 0.1 g cm~^, typical of diffuse star-forming 
clouds such as Perseus and Ophiuchus. The second has 
S = 1.0 g cm~^, typical of regions of massive star for- 
mation in the Galaxy. The third has S = 10 g cm~^, 
an extremely high surface density found only in clusters 
near the Galactic center and in extragalactic super star 
clusters. However, this type of star-forming environment 
is believed to have been more common earlier in cosmo- 
logical evolution. We call the runs L, M, and H, for low, 
medium, and high column density. We summarize the 
setup of the three runs in Table m We run each simula- 
tion for a time t = O.Gig , where 



37r 
32Gp 



(17) 



is the free-fall time computed at the mean density p = 
iM / {AttR^) of the cloud. By this point the differences 
between the runs are clearly established, and the fraction 
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TABLE 1 

Simulation Parameters 



Name 


M (Mq) 


S (g 


R (pc) 


cTy (km s -'-) 


iff (kyr) 


ibox (pc) 


L 




Axo (AU) 


Nl 


AxL (AU) 


L 


100 


0.1 


0.258 


1.29 


217 


1.95 


8 


256 


1573 


65,536 


6.14 


M 


100 


1.0 


0.081 


2.30 


38.6 


0.489 


6 


256 


393 


16,384 


6.14 


H 


100 


10.0 


0.026 


4.08 


6.86 


0.160 


5 


168 


197 


5,376 


6.14 



Note. — Col. 7: linear size of computational domain. Col. 8: maximum refinement level. Col. 9: number of cells per linear 
dimension on the coarsest level. Col. 10; linear cell size on the coarsest level. Col. 11-12: same as col. 9-10, but for the finest 
level. 



of the collapsed mass in the most massive star asymp- 
totes to a constant value (see below). 

We emphasize that we have chosen the numerical setup 
so that the runs are, as much as possible, simply rescaled 
versions of one another. The initial conditions have iden- 
tical density structures, virial ratios, and velocity fields, 
and the refinement criteria and peak resolution is the 
same in every run. We simulate each cloud for the same 
number of free-fall times. The homology between the 
runs is broken only by the infiuence of radiation. Radia- 
tion fixes the gas temperature (and thus causes the initial 
Mach numbers of the runs to vary slightly), and, much 
more importantly, the very different optical depths of the 
different clouds cause them to respond differently once 
stellar feedback begins. Thus we expect any difference 
between the three runs to be almost entirely dictated by 
their differing response to stellar radiative forcing. 

3. RESULTS 

3.1. Collapse Morphology 

We show the large-scale evolution of runs L, M, and H 
in Figurejl] As the plot shows, and as expected, the three 
runs are essentially homologous on large scales. Once 
we account for the scaling of cloud radius and surface 
density between the runs, the only noticeable difference is 
a slight trend toward increasingly filamentary structures 
as E increases. This is easily understood as a radiative 
effect. As noted above, runs with higher S require larger 
velocity dispersions to maintain initial virial balance, 
while the sound speed is fixed by radiative effects. Thus 
the Mach number of the initial turbulence increases with 
E, and this produces the observed increase in filamentary 
structure. 

In Figure [2j we show the simulations at the same times 
as in Figure m but now zoomed in to a small region 
centered on the most massive star present, or on the 
origin before stars form. In the left panel of Figure [2] 
the images are scaled homologously, so that the regions 
shown all have a length equal to 10% of the initial cloud 
radius, and the color scale is in units of surface density 
divided by initial mean surface density. In the right panel 
the scaling is physical, so that all the plots show a region 
of fixed physical size, using a column density scale in 
fixed rather than normalized units. 

In contrast to the large-scale homology seen in Fig- 
ure [l] on small scales the runs rapidly diverge. Devia- 
tions from homology start to appear around 0.2ift, and 
by OAtff any visual similarity between the runs is com- 
pletely gone. Progressing from low to high S, runs are 
characterized by increasing surface densities even when 
normalized to the initial surface density. In run H the 
predominant structure is a single large disk concentrated 
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Fig. 1. — Column density in simulations L, M, and H (left to 
right column) at times running from t = to t = 0.6tfi {top to 
bottom row). The color scale is normalized to the initial mean 
column density Eq = 0.1, 1, and 10 g cm~^ for runs L, M, and H, 
respectively. 
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around a single central object. In run M we have a mas- 
sive binary with two circumstellar disks and a larger cir- 
cumbinary disk. Finally, in run L the disks are much 
smaller and less dense, and they have mostly depleted 
by the final time shown. 

3.2. Fragmentation and Star Formation 

The higher E runs also fragment less, producing fewer, 
more massive stars than the low E runs. We show this 
in Table [5] and Figure [3] The total mass in stars M*_tot 
at any given time (normalized to tg) is very similar from 
run to run, varying by less than 10% from run L to run H 
at times > 0.2tff . At the end of the simulation, each run 
has a total of 15 Mq of stars. This reflects that the star 
formation rate in the simulations is driven by large-scale 
flows that are changed very little by radiation feedback. 
However, the pattern of fragmentation is quite different. 
The mass of the most massive object M^^max in the three 
runs begins to diverge at 0.2 — 0.3tg, and thereafter it 
increases much more rapidly in run H than in run L. At 
the final time the difference in mass is nearly a factor of 3, 
with run L having a maximum stellar mass below 6 Mq , 
and run H reaching 15 Mq. In all the runs the fraction 
of the total stellar mass in the most massive object /max 
asymptotes to a roughly constant value after O.Btff . The 
asymptotic value ranges from /max 0.35 in run L to 
/max ~ 0.9 in run H. In effect the initial gas cloud in run 
H is like a single massive protostellar core that forms 
one massive star plus a few small secondaries, while the 
cloud in run L instead forms a small cluster that does 
not include any massive stars. Run M is intermediate. 

The difference in runs is even more apparent if we focus 
on a single time. Figure [4] shows the cumulative distri- 
bution function of stellar mass at t = O.Gts in each of the 
runs. In run L the distribution function rises relatively 
smoothly between 1 and 5 Mq, so the system consists 
of a number of stars of roughly comparable mass. In 
contrast, in run M 90% of the mass is in just two stars 
th at form a binary system, a result very similar to that 
in |KKM07| which used an initial surface density S = 0.7 
g cm~^. In run H a comparable fraction of the mass is 
in a single star. Since the system in run L involves a 
number of stars of comparable mass, these are unlikely 
to wind up as a bound star system. The system in run 
M, on the other hand, is stable and is likely to remain 
a bound binary. Thus in both runs M and H, the result 
is that most of the mass goes into a single star system, 
while in run L the mass will end up divided into several 
star systems. 

3.3. Thermal Structure 

The difference in morphology, fragmentation, and star 
formation is easy to understand if we examine the tem- 
perature structure of the gas. In Figure [5] we show 
the column-density weighted temperature over the same 
small-scale regions as in Figure |2] Clearly the tempera- 
ture distribution in the gas in the runs is even less ho- 
mologous than the density structure. Moreover, the dif- 
ferences in temperature begin to appear at earlier times. 
At t = O.ltff the density plots are nearly indistinguish- 
able once the homologous scaling is removed (left panel 
of Figure [2| , while the differences in temperature are al- 
ready obvious. It is important to point out that at this 



point there no stars present that are producing signifi- 
cant amounts of power via nuclear luminosity. The most 
massive star present in any of the runs is only 0.37 Mq, 
and its luminosity is entirely driven by accretion. The 
difference in temperature is therefore solely due to the 
two factors pointed out by |KM08[ compared to run L, 
run H has both a higher density to produce higher ac- 
cretion rates and thus higher accretion luminosities, and 
a higher optical depth to more effectively trap the radi- 
ation that is produced. 

We can understand the difference in the thermal struc- 
ture of the runs more quantitatively, and explore how this 
difference is likely to infiuence fragmentation, by exam- 
ining the relationship between temperature and density 
in each of the runs. Figure [6] shows the locus occupied by 
the clouds in runs L, M, and H at different times in the 
simulation in the plane of log density versus log temper- 
ature. The color represents the mass density at a given 
point, and the contours indicate, from lowest to highest, 
regions in the plane containing 99.9%, 99%, 90%, and 
50% of the gas mass in the computational domain. 

The figure demonstrates how different the thermal 
structure of the gas is in each of the runs. In run L, the 
great majority of the gas is near the background tem- 
perature of 20 K at all times. Even at the final time 
only 10% of the gas mass is at noticeably elevated tem- 
peratures, and less than 1% of the mass is heated above 
100 K. In contrast, in run H the heating is much more 
extensive. Even at t = O.lift, when the only source of 
heating is the accretion luminosity of a 0.37 Mq star, 
the contours containing 50% and 90% of the mass are 
noticeably elevated above the T = 20 K line. Deuterium 
burning in the star begins shortly before 0.2tff, and by 
the final time deuterium burning and Kelvin-Helmholtz 
contraction (the star has not yet reached the main se- 
quence) provide enough luminosity to keep all the mass 
at temperatures > 50 K. Run M is intermediate, with 
small but significant fractions of the cloud mass reach- 
ing elevated temperatures at early times, and more mass 
becoming heated as the run progresses and the stars be- 
come more luminous. This is consistent with the analytic 
predictions of KM08, who find that a column density of 
1 g cm^^ constitutes the rough line between clouds that 
do and do not experience significantly elevated temper- 
atures over much of their mass as a result of trapped 
accretion luminosity. 

It is important to point out that the temperature does 
not need to rise to the point where the Jeans mass is 
above 100 Mq in order to inhibit fragmentation. Indeed, 
such a rise in temperature would be sufficient to halt col- 
lapse of the core entirely. Instead, the heating prevents 
fragmentation by creating an environment where the ef- 
fective equation of state with 7 = 1-1- dlogT / d\og p > 1 
throughout the bulk of the cloud mass. Examining Fig- 
ure |6j we see that the region containing 90% of the cloud 
mass (the third contour from the outermost one) is al- 
most perfectly horizontal in run L at all times, so 7 « 1. 
In runs M and H, on the other hand, this region has a 
slope ^ 0.2 — 0.3 in the log p — log T plane at al l times of 
0.2tff or more, similar to the result obtained by |KKM07| 
indicating that t he effect ive equation of state is closer to 
7 = 1.2 — 1 .3. As Larson (2005) points out, and the sim- 



ulations of Ja ppseiT et al.| ( j2005p confirm, fragmentation 
is likely as long as the effective equation of state for the 
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Fig. 2. — Column density in simulations L, M, and H {left to right column) at times running from t = to t = 0.6% {top to bottom 
row). Symbols indicate stars, with the type of symbol indicating the stellar mass. Low mass stars (M« = 0.05 — 1 Mq) are indicated by + 
signs, intermediate mass stars (M* = 1 — 8 Mq) by X signs, and massive stars (M* > 8 Mq) by filled circles. Left: the region shown is a 
0.1-R X 0.1-R box centered on the most massive star, or the origin if no stars are present, and the color scale is normalized to the initial mean 
column density S, where R and E have the values given in Table[l]for runs L, M, and H. Right: the region shown is a 3000 AU X 3000 AU 
box centered on the same point as in the left panel, and the color scale is in the same physical units for every run. 



gas is 7 < 1, and is unlikely when 7 > 1. In runs M and 
H, there is no significant gas mass with 7 < 1, which is 
why fragmentation is suppressed. 

Finally, we emphasize that these phenomena cannot 
be correctly captured by analytic equations of state that 
are based on either a baroptropic or and optically thin 
cooling as sumption. Exa mpl es of suc h equations of state 
from Dobbs et al. ( 2005 1 and Larson ( ,2005 ) are shown in 
Figure |Sj and they clearly do not even come close to re- 



producin g the result s with radiative tran sfer, a point also 



made by Boss et al.'(2000), Krumholz (2006|, Krumholz 



et aL| ( |20U7ap , and |Ottncr et al.^ (2009bj . "Sny such ap 
proximation would give the same temperature-density re- 
lation for all three of our simulated clouds, while clearly 
the results are different at different times and for differ- 
ent initial cloud column densities. [Urban et al. (20101 



reach the same conclusion for lower- density, larger-scale 
clouds based on their simulations. 
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TABLE 2 
Stellar Content versus Time 



i/*ff ^* Mf^tot Af,,niax /max M,_tot Af,,max /max -^V, A/»_tot A/.^max /max 
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0.18 


1.00 
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0.29 


1.00 
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0.37 


0.37 


1.00 


0.2 
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0.95 


0.86 


0.90 


3 


1.55 


1.16 


0.75 
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1.91 


1.86 


0.97 


0.3 


3 


5.06 


2.54 


0.50 
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5.73 


2.56 


0.45 


3 


5.60 


5.13 


0.92 


0.4 


7 


8.10 


2.93 


0.36 


7 


9.61 


5.07 


0.53 
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8.16 


7.10 


0.87 


0.5 


7 


11.55 


4.06 


0.35 
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12.44 


7.18 


0.58 


2 


10.92 


10.83 


0.99 


0.6 


8 


15.58 


5.75 


0.37 


7 


16.42 


8.77 


0.53 


9 


14.96 


13.41 


0.90 



Note. — Col. 1: run time. Col. 2: number of stars present in run L. Col. 3: total mass of stars in run L. Col. 4: mass of 
largest star in run L. Col. 5: fraction of total stellar mass in largest star. Col. 6-9 and 10-13: same as columns 2-5, but for runs 
M and H. 

capture the difference in optical depth between runs, but 
it would not capture the effect that higher density runs 
produce higher accretion rates and thus higher accretion 
luminosities from the embedded protostars. The analytic 
models of KM08] suggest that both effects are of compa- 
rable importance. 

4. DISCUSSION 

4.1. The Massive Star Fraction 

Our results demonstrate that the amount of fragmenta- 
tion that a cloud undergoes is likely to depend strongly 
on its surface density, and this has important implica- 
tions for where we expect massive stars to form. Con- 
sider a protostellar core, an object with a mass of a few 
tenths to a few hundreds of Mq that collapses to make 
one or more stars. L ow mass cores do not have significant 
internal turbulence ([ An dre et al.||2007[ |Kirk et al^|2007 




0.5 0.6 



Fig. 3. — Star formation histories in runs L, M, and H. Top: 
total stellar mass Af«^tot {thick lines) and mass of most massive 
star Aft^max {thin lines) versus time, as indicated. Bottom: frac- 
tion /max of total stellar mass in most massive star versus time. 
Sharp jumps represent mergers between a central massive star and 
a smaller star. 
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grounds ( Uttner etaf. 2008b) . Consequently, while they 
may fragment into a binary like run M, we expect most 
of the stellar mass they produce to end up in a single 
star system. The overall efficiency of turning gas mass 
into stellar mass is expected to be e « 1/3 rather than 
e = 1 as a result of mass e.]ection by protostellar outflows 

Alves et al.||2007| [Enoch et aL 



(Matzner & McKee[[2000 



2DDSU 



Fig. 4. — Fraction /(< M) of total stellar mass contained in 
stars with mass < M as a function of M, for each of the runs at 
time t = 0.6tff. 



Although we have not tested the Bate (2009') approach 
of omitting radiation from stars and including only ra- 
diative emission by gas on size scales resolved by the 
computation (which are much larger than stellar scales 
in both Bate's calculation and ours), it seems unlikely 
that this approximation could succeed either in the case 
of clouds with differing initial column densities. It would 



Tn contrast, as pointed out by McKee fc Tan[ (|2003), 
massive cores such as those we simulate are turbulent, 
since their masses are many times the thermal Jeans 
mass. We find that, at low S, such cores will fragment so 
that their stellar mass is divided among many star sys- 
tems. Figure [7| gives a more detailed picture of how this 
happens in run L. There are five fragments whose masses 
appear to be asymptotically approaching fixed fractions 
of the total stellar mass, ranging from 11% to 37%, plus 
three more smaller stars whose masses appear to have 
reached nearly fixed maxima, and that are therefore de- 
clining with time in the total fraction of stellar mass that 
they represent. 

We can use this result to make a toy model for how 
the fraction of the stellar mass that is in massive stars 
is likely to vary with surface density. We begin from the 
observation that stars form from cores that have a mass 
distribution with the same functional form as the IMF, 
so that the IMF is set at the phase w hen gas fragments 



into protostellar cores ( Motte et al.|19 98 ; Tcsti & Sargent 
1998; Johnston e et al.||2001[ 'Unishi ct al. 200 2i |Beuther 
et al. 2004: Reid & Wilson|,2005, 2006a b, Alves et al 
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Fig. 5. — Same as Figure [2] except that the plots show column density-weighted temperature, defined as f pTdz/ f pdz. The color 
scales are the same in both the left and right sides. 



2007 



2008 ^ 
model' 



Nutter fc Ward-Thom psonI 1 2007] jSimpson et al 



Enoch et al. 2008| Kathborne et al. 
this core mass tunction (CMF 



2009[). We 



Chabrier 



usmg a 

(120051) stellar system IMF shifted to higher mass by a 
factor of 3 to account for the mass that is ejected by 
outflows: 



dUc 



din nic 
with nTc 



(In rric — ln nic) 
-1.3 



Aexp 
B ( -^^) 

\mbioak J 

0.75Mq, mbreak = SM© , and a — 0.55 



m,c < mbreak 

rric > mbreak, 
(18) 



Our values of rUc and mbreak are chosen so that, when 
e = 1/3, the stellar IMF will have a peak at 0.25 Mq 
and will break from log normal to powerlaw form at 1.0 
Mq , in agreement with |Chabrier| s best fit to observa- 
tions. The normalization factors are related by i? = 
Aexp[— (luTObrcak — In m^)^/(2(T^)] . Theoretical models 
are able to explain this distribution of core masses as 
arisin g naturally from th e proper ties of supersonic turbu- 
lence ( Padoan & Nordlund 2002 ; Hennebelle & Chabrier 
2008 1. In this picture, the final stellar IMl'' is simply 
the convolution of the core mass function (CMF) with 
a simple function that maps core mass to stellar mass. 
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Fig. 6. — Temperature-density relation for runs L, M, and H 
(left to right columns) at the same times as in Figure IT] except 
that we omit time t = 0.0. Colors show the mass density per 
square dex in the log p — log T plane. Contours from outermost to 
innermost indicate the region in the log p — log T plane containing 
99.9%, 99%, 90%, and 50% of the total gas mass at that time (i.e. 
not including the mass in stars). In some runs only the outermost 
contours are visible because the inner ones form a thin line near 
T = 20 K. The linear feature extending to high temperature and 
low density represents cells that contain a mix of cloud and hot 
ambient medium; such cells never contain a significant fraction of 
the cloud mass, as indicated by the contours. The solid line is 
the barotropic curve of iDobbs et al.| ( |2005l), which has a constant 
temperature T = 20 K at low densit y. I'he dashed line is the 
optically thin cooling approximation of |Larson| | [2005| ). 



and w hich must be nearly mass-independent. [Clark et aL| 
( 2007 1 suggest that this correspondence will be disrupted 
if the core free-fall time is mass-dependent. However, ob- 
served cores do no t have mass-dependent free-fall times 
(Andre et al."20G7), and theoretical models predict that, 



contrary to Clark et al.|' s assumption, co re free-fall time 
depends on mass at most very we akly (McKee & Tan 
2003| [Hennebelle fc Chabrier]|2009[ ) as a result of turbu- 
lent support. 

To make a quantitative model of how fragmentation in 
low S regions will affect the IMF, we consider two ex- 




FlG. 7. — Mass (top panel) and fraction of the total stellar mass 
(bottom panel) in all stars as a function of time in run L. Each line 
represents an individual star. 

treme scenarios that bracket the outcome of run L. The 
first is that, in regions of low S, cores with mass ttLc above 
some minimum fragmentation mass Wfrag produce only 
a single massive star with a mass m* = mc/9, i.e. 2/3 of 
the core mass is ejected by outflows, 1/3 of what remains 
goes into the largest star, and the remaining 2/3 goes into 
low mass stars with m* < mfrag/3. The alternative pos- 
sibility is that the cores with mass TOc > Wfiag fragment 
into rtfrag stars of equal mass = mc/(3nfi.ag), with 
the factor of 3 to account for ejection by outflows. Based 
on the outcome in run L, we adopt nfrag = 5 for our 
toy model, although of course in reality we expect that 
the number of fragments and their mass distribution will 
vary stochastically. 

In the first scenario, if the total mass of cores with 
masses between and mc+dnic is given by dnd dhmric, 
the corresponding mass of stars with masses between 
— rad'i and m* + dm* is given by dn^,/d\Tim^ — 
( 1/9) dric/fi In TOc for any mass > mftag- In the sec- 
ond scenario, we instead have to* — rric/ (3nfrag) and 
dn^/dlnm^ — (l/3)dnc/d\nmc, since all of the core 
mass that is not ejected by outflows (i.e. 1/3 of it) 
goes into stars of mass TOc/(3nti.ag). Finally, in high S 
regions where massive cores do not fragment, we also 
have rfn^/dhiTO* = (l/3)c?nc/dlnTOc, but now the stel- 
lar mass is related to the core mass by m* = TOc/3 rather 
than TO* = TOc/(3nfrag). 

Given these relations, we can compute the fraction of 
all stellar mass that is contained in stars with masses 
greater than to*, which we denote /(> m*). We adopt 
a minimum stellar mass TO * niin ~ 0.01 and a max- 
imum TO* max = 120 Mq (Figer 20051. First consider 
regions of high S where massive cores do not fragment 
and stellar and core masses are related by to* = TOc/3. 
In such a region, we have 



/(> TO*) = 



J3. 



d In mc 



dnir 



3m* 
3m* 



din 



drrir 



fH{>m,). (19) 



This equation simply states that the fraction of stellar 
mass in stars larger than to* is the same as the fraction 
of core mass in cores larger than 3to* . For low E regions, 
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Fig. 8. — Top panel: fraction of mass /(> m») in stars with 
mass greater than m*, for regions of high surface density (thick 
Une, fu, Equation |19| and regions of low surface density (thin 
hnes, fL,i, /l,2, fh.ia, and fL,2a, corresponding to Equations [20| 
— 123| . Bottom panel: ratio of /l(> fHi> ti*) massive star 

mass fraction in low S regions to that in high S regions. As in the 
top panel, Jl = /l,2, fL,la, or fL,2a, as indicated. 

in the first scenario we instead have 



/(> ™*) 



J9m* 
■J 3m*, n 



d In rrif. 

dric. 
d In nic 



drric 
drrir 



- /l,i(> ^>^) 



(20) 

for stellar masses to* > TOfrag/3. Note the factors of 1/9 
in the numerator and 1/3 in the denominator, reflecting 
that, in this scenario, only 1/9 of the mass in a core of 
mass nic is incorporated into a star of mass m* = rric/Q, 
but that 1/3 of the core mass goes into stars overall. If 
we instead adopt the second scenario, where 1/3 of the 
core mass goes into Ufrag stars of mass to* = TOc/(3ntrag), 
we obtain 



/(> m*) 



3n 



frag? 



d\n 



drrir 



3m* ^ 



d In nic 



drrir 



/l,2(> ™* 



(21) 

A final complication is that, in both scenarios, we have 
implicitly assumed that the core mass function extends 
to infinity, or at least to 3?T.fragTO*^iiiax = 1800 Mq. This 
is not necessarily the case - the origin of the observed 
cutoff in the IMF is not understood, and one possible 
explanation for it is that there simply are no protostellar 
cores whose mass is larger than a few hundred Mq that 
are capable of collapsing to single stars even in regions 
of high surface density. As a simple example of how this 
would change the results, suppose that there is a maxi- 
mum core mass m^max = 3TO«_i„ax — 360 Mq, sufhcient 
to make a 120 Mq star if the core does not fragment. In 
this case the mass fractions and /l_2 we have just 
computed are modified to 



/L,la(> m^*) - 



/L,2a(> m^) - 



pc,„.x /IN 
J9m, \9l 



dnc 



i In nic 



drrir 



rnic 
J 3m 



3nfra 



£ ( 1\ dn 
„ V3/ din 7 

= n \ _d, 
^, V3/ din 



drric 
- drrir 



(22) 



(23) 



Evaluating the functions fn, /l.2, /L,ia, and 

/i,2a gives the results shown in Figure [Sj As the plot 



shows, reduced star formation efhciency due to the frag- 
mentation in massive cores of the sort we have found 
can reduce the fraction of total stellar mass in high mass 
stars by a significant amount. The minimum reduction in 
massive star fraction, by 40%, is for /2.Lj corresponding 
to the scenario where massive cores fragment into a few 
equal mass objects and the core mass function extends to 
infinity. In this case the full mass of cores that fragment 
is still available to make massive stars, and the massive 
star fraction declines only because stars of mass to* must 
be produced by cores of mass 3nfragTO* in regions of low 
E rather than by cores of mass 3to* in regions of high S, 
and the total mass of cores available is lower for higher 

TOc. 

Any other scenario gives a much more significant re- 
duction in the massive star fraction in regions of low sur- 
face density. In scenario where fragmenting cores 
produce one massive star of mass TOc/9 and only low 
mass stars otherwise, the reduction is by 76%. The re- 
duction is partly for the same reason as for /l,2, and 
partly because all the mass in the massive core that does 
not go into massive stars still goes into low mass stars. 
Finally, if the core mass function has a cutoff, the reduc- 
tion in massive star fraction is even more dramatic. In 
this scenario fragmentation means that a core of mass 
TOc = 9to* (for fh.ia) or of mass TOc = 3nfragTO* (for 
fL,2a) is required to make a star of mass to*, and if this 
exceeds the maximum core mass then no stars of mass 
TO* can form in regions of low E. 

Regardless of which scenario is ultimately correct, our 
results show that in regions of low surface density we ex- 
pect a significant decline in the stellar mass fraction in 
massive stars compared to a canonical IMF. The reduc- 
tion is anywhere from a factor of 1.7 in the most con- 
servative scenario to a very large factor in more liberal 
scenarios. This provides numerical confirmation of the 
hypothesis advanced by |KM08| that there is an effective 
threshold for massive star formation. 

4.2. Implications of Environmental Variation in the 

IMF 

A variable IMF has numerous implications on scales 
ranging from the sub-galactic to the cosmological, and a 
full exploration of them is beyond the scope of this paper. 
Moreover, a full exploration of this topic would require 
a larger parameter study than the one we present here, 
allowing a full exploration of how fragmentation and the 
stellar mass function vary with environment. Nonethe- 
less, we can identify some important implications. 

On the scales of star clusters, preferential formation 
of massive stars in regions of high surface density sug- 
gests that clusters are likely born mass-segregated, with 
more massive stars forming in the center where the sur- 
face density is highest. Outer parts of the cluster, where 
the surface density drops below ^ 1 g cm~^, should form 
preferentially low mass stars. There is some evidence 



for such pr imordial mass segregation in Orion (Huff & 
|Stahler°2006), but debate continues about whether mass 
segre gation in clusters in gene ral is a result of forma- 



tion ( [Bonnell fc Davies|[T998 
during~Ehe hrst tew IVLyr of c. 



or dynamical 



2006 
thoug 



uster lifetime (Tan et al 



processes 



McMillan et al. 2007). Both may occ ur to gether 



r_ recent observational fFiiresz et al."2008! 'Tobinl 



et al.pOOgj ) and theoretical ( jOftncr 



et al.„2008b, 2009al 
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work pointing out that stars are born sub-virial with re- 
spect to their parent clouds (even if the clouds themselves 
are virialized) suggests that dynamical segregati on may 



be much st ronger than earlier estimates found (Allison 
et al |2009[ ). In very dense clusters, primordial mass seg- 



regatlon may also strong ly influence the dynam ical evo 
lution of the cluster (e.g. Chatterjee et al.||2009 ). 

Conversely, we do not expect a correlation between the 
IMF and cluster mass beyond the trivial one expected 
from finite sampling of a universal IMF, since there does 
not appear to be a strong correlation between pr oto clus- 
ter g as cloud masses and surface densities (iFall et al. 



20101. A correlation between cluster mass and the IMl*' 



has has been su ggested by some models (e.g. Weidner 
fc Kroupa||2006 ), based on idea that massive stars form 



via a process m which all gas clouds fragment down to 
the small initial Jeans mass, but that some subseq uently 



grow through c ompetitive Bondi-Hoyle accretion (Bon 



nell et al. 2004). However, the premise on which these 
models are based - fragmentation of all clouds down to 
the Jeans mass at the initial, low cloud temperature - 
clearly fails when radiation is included. While the pro- 
cess that goes on in run L might roughly be described as 
competition, there is clearly no competition in runs M 
or H, where radiation ensures that most of the proposed 
competitors never form in the first place. 

Observations remain divided on whether there is a 
correlation between cluster mass and maximu m stellar 

and Weidner 



mass. 



2006h 



Weidner fc Kroupa| ( |2004, 
et all (2009) claiii i to d etect one, while luey et al. ( 2004^ , 
Ermegrccnf( |2006[ ), and [Parker fc Goo dwin (2007) argue 



tFiat there no correlation is present, de Wit et al.| (2004 



2005 1 find that 4± 2% of galactic O stars formed outside 



of a cluster of significant mass, which is consistent with 
the models presented here (for example runs M and H 
form effectively isolated massive single stars or binaries) , 
but not with the proposed cluster-stellar mass correla- 
tion. 

No comparable studies have been conducted to search 
for systematic variation of the IMF with surface den- 
sity, which we do predict. Such studies are likely to be 
even more challenging than those searching for a correla- 
tion with cluster mass, because cluster surface densities 
evolve very rapidly once star formation ends and the re- 
maining gas is expelled. For this reason any search for a 
correlation between IMF and surface density would have 
to target clusters that are still embedded in their parent 
gas clouds, for which the cloud surface density can be 
measured. Unfortunately this renders optical observa- 
tions, the most common method for determining stellar 
masses, impossible, since a surface density E = 1 g cm~^ 
corresponds to Ay ~ 200. Instead other methods of es- 
timating populations of low and high ma ss stars, such 



as x-ray obse rvations, would be required (Krumholz & 



McKee 



2008). 



On larger scales, our ability to make predictions is lim- 
ited by our lack of a theoretical model capable of connect- 
ing the surface densities measured over the small scales 
of star-forming regions where fragmentation occurs (~ 1 
pc) to those averaged over much larger areas of galac- 
tic disks (~ 1 kpc). Regions of high surface density 
where massive stars form are clearly found preferentially 
in regions of high galactic surface density such as spiral 
arms and near galactic centers, while low surface density 



clouds such as Taurus are found in regions of lower sur- 
face density. However, there is no clear one-to-one map- 
ping from large to small scales. Nonetheless, it seems 
clear that our results do predict a suppression in the for- 
mation of high mass stars in regions where the galactic 
surface density is low, such as dwarf galaxies and the 
outer parts of spiral galac tic disks . Such a correlation 



may have been observed (Boissier et al. 2007 Meurer 



et al.|2 009 Lee ct al. 2009), although it too remains con 



troversial ( )Boselh et al.||2009| . 



5. SUMMARY 

We report the results of a series of adaptive mesh re- 
finement radiation-hydrodynamic simulations of the col- 
lapse of a massive gas cloud using our code ORION. In 
these simulations, the initial density and velocity struc- 
ture, the initial virial ratio, and all numerical aspects 
of the runs are held constant, but the clouds are scaled 
to different initial surface densities, leading them to re- 
spond differently to the radiation feedback produced by 
the stars that form in the cloud. We find that this dif- 
fering response leads to dramatic differences in the ways 
the clouds fragment. A cloud with an initial surface den- 
sity E = 0.1 g cm~^, typical of low-mass star-forming 
regions such as Taurus or Perseus, fragments strongly, 
such that it produces a number of stars of comparable 
mass and puts only a small fraction of its initial mass into 
the most massive star. No massive stars form in this run. 
In contrast, runs with initial surface densities of 1 and 10 
g cm~^, typical of Galactic massive star- forming regions 
and extra-galactic super-star clusters, respectively, frag- 
ment much less. The run with E = 1 g cm~^ puts most 
of its mass into a single massive binary system, while the 
one with E = 10 g cm~^ ends with 90% of the stellar 
mass in a single star. 

We show that these differing outcomes can be under- 
stood in terms of the way that radiation feedback from 
the stars forming in the cloud affect its subsequent frag- 
mentation. The higher surface density clouds are char- 
acterized by higher accretion rates, leading to higher ac- 
cretion luminosities at early times. Furthermore, their 
higher optical depths trap the resulting radiation more 
effectively. The net effect is that radiation feedback raises 
the temperature and thus the Jeans mass over a signif- 
icant fraction of the cloud mass in the highest surface 
density runs, while affecting only much smaller regions 
when the surface density is low. This leads to an increas- 
ing suppression of fragmentation as the initial surface 
density rises. 

Our results suggest that the stellar IMF need not be 
universal between regions of low surface density (E <C 1 
g cm~^) and those of high surface density (E > 1 g 
cm~^). In the former, even if turbulence creates the same 
mass spectrum of initial protostellar cores as in the lat- 
ter, these cores will fragment during collapse, producing 
small clusters rather than individual star systems. This 
effect can dramatically reduce the fraction of the mass of 
a stellar population in stars with masses > 10 Mq. 
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